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Nomenclature 


u 

=  Solution  vector 

E,  F 

=  Fluxes 

u,  V 

=  Cartesian  velocity  components 

P 

=  Density 

P 

=  Pressure 

e 

=  Internal  energy 

E, 

=  Total  energy 

ttr,  TXy%  Tyy 

=  Shear  stress  components 

q*.  qy 

=  Heat  Flux 

Re 

=  Reynolds  number 

Pr 

=  Prandtl  number 

M 

=  Mach  number 

T 

=  Temperature 

(X,y) 

=  Cartesian  coordinates 

(V.  <p) 

=  Streamfunction-potential  function  coordinates 

J 

=  Jacobian 

(xv,  xv...) 

=  Grid  metrics 

(x„  yj 

=  Grid  speed 

Sv 

=  location  of  the  first  grid  line  above  the  wall 

=  Strength  of  vortices 

T 

=  Clustering  parameter  in  <p  direction 

P.Q 

=  Clustering  parameters  in  \\i  direction 

9 

=  Angle  between  grid  lines 

S 

=  Metric  tensor 

St 

=  Strouhal  number 

I.  Introduction 

Von  Mises  (Ref.  1)  was  the  first  to  introduce  streamfunction  as  a  coordinate  (SFC).  In  his  study,  he  used 
the  streamfunction  as  the  coordinate  across  boundary  layer  while  the  transverse  Cartesian  coordinate  was 
unchanged.  A  limitation  of  this  approach  is  the  degeneracy  of  the  grid2  near  the  leading  edge  of  the 
airfoil.  A  number  of  modifications  have  been  proposed  to  deal  with  this  issue2'5.  The  advantage  of  the 
SFC  methodology  in  CFD  lies  in  the  fact  that  the  streamfunction  serves  both  the  purpose  of  grid 
generation  and  description  of  motion  for  steady  flows.  The  SFC  approach  has  found  applications  in  a 
wide  variety  of  fields  including  two-dimensional  channel  design6,7,  nozzle  inverse  flow  design 
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problem8, inverse  design  of  aerodynamic  shapes9,10  and  turbomachinery  flows,  to  name  a  few.  In  the 
aerodynamics  of  moving  airfoils,  the  generated  SFC  grid  can  move  with  the  airfoil  and  keep  the  boundary 
layer  well-resolved.  This  feature  is  crucial  for  accurate  modeling  of  flow  separation. 

SFC  allows  for  the  use  of  economic  single-grid  high-order  finite-difference  type  schemes  for  complex 
geometry.  In  terms  of  computational  time  per  grid  point,  this  method  is  advantageous  compare  to  multiple 
grids  or  unstructured  grids  combined  with  high-order  schemes  and  interface  matching  conditions. 
Although  the  SFC  is  widely  used  in  combination  with  lower-order  numerical  schemes  in  CFD,  the 
straightforward  implementation  of  SFC  create  rather  coarse  grid  at  the  vicinity  of  stagnation  points  that 
smears  high-order  numerical  computations  of  aeroacoustics  and  unsteady  aerodynamics  of  flapping 
airfoils.  The  application  of  the  SFC  methodology  for  numerical  computations  with  high-order  of 
approximation  is  considered  in  the  present  report. 

Though  the  main  application  of  SFC  was  in  the  field  of  inverse  flow  device  design,  it  has  found 
application  in  other  fields  too.  Finnigan  (Ref.  1 1-13)  used  the  SFC  approach  to  study  turbulent  shear  flows 
and  chaotic  systems.  The  method  has  also  found  application  in  the  field  of  two-phase  flows14  and  flow 
through  porous  medial5.Chakravarthy16  used  the  streamlines  and  potential  lines  of  an  incompressible 
solution  to  generate  a  grid  for  transonic  nozzle  flow.  Adamczyk17  generated  O-type,  C-type  and  split 
conformal  grids  for  turbomachinery  flows  using  an  integral  formulation  similar  to  the  Symm 
construction18  to  generate  orthogonal  grids. 

The  streamfunctions  are  defined  only  for  two-dimensional  flows;  nevertheless,  applications  can  be 
found  in  literature  where  the  SFC  has  been  applied  to  3-D  problems.  Anderson19  generated  a  three- 
dimensional  coordinate  system  by  rotating  an  axisymmetric  2-D  duct  about  its  axis.  Briley20  produced  the 
grid  for  a  3-D  duct  by  stacking  2-D  systems.  Hounjet 21, 22  generated  2-D  and  3-D  grids  using  the  panel 
method  solution  of  the  velocity  potential.  The  method  was  used  to  generate  O-type  and  multi-block  OH- 
type  grids  about  transverse  cross-sections  of  transport  type  aircrafts.  As  opposed  to  O-grids  and  C-grids, 
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the  SFC  coincides  with  Cartesian  coordinates  in  the  far-field  that  makes  straightforward  the  use  of 
outflow  boundary  conditions. 

The  problem  with  the  SFC  is  that  a  uniform  distribution  of  grid  points  in  the  computational  domain 
produces  coarse  cells  close  to  the  stagnation  point.  For  the  only  velocity  scale  in  the  computational  set-up, 
the  velocity  is  slow  at  the  vicinity  of  stagnation  point  and  the  coarse  grid  does  not  introduce  significant 
error.  However,  this  problem  is  critical  for  unsteady  aerodynamics  including  flapping  airfoils  because  of 
the  presence  of  another  velocity  scale  such  as  speed  of  moving  airfoil.  Indeed  the  speed  of  moving  airfoil 
is  large  for  higher  frequency  of  airfoil  pitching  and  plunging.  The  error  introduced  by  coarse  grid  may 
harm  the  effort  to  apply  high-order  schemes  to  problems  of  flapping  airfoils. 

To  circumvent  this  issue,  the  following  combination  of  algorithmic  steps  has  been  implemented  in  the 
current  study.  A  hyperbolic  stretching  function23  was  used  for  grid  refinement  close  to  the  solid 
boundaries.  The  degree  of  stretching  will  be  chosen  carefully  to  make  sure  that  the  grid  does  not  become 
extremely  coarse  in  the  far  field.  The  spatial  differencing  is  performed  using  a  non  uniform  compact 
scheme24,  25.  To  remove  the  oscillations  due  to  the  grid  non-uniformity,  an  implicit  Pade-type  filter26  is 
employed  in  the  current  study.  The  explicit  time  marching  is  performed  using  a  Runga-Kutta  algorithm, 
implemented  in  the  low-storage  form27.  A  simple  energy  transfer  and  annihilation  method28  is  used  at  the 
domain  boundaries  for  absorbing  the  outgoing  energy. 

The  governing  equations  are  listed  in  Sec.2.  Sec.3  discusses  the  SFC  method.  The  numerical  algorithm 
is  discussed  in  Sec.  4.Numerical  results  are  discussed  in  Sec. 5.  In  Sec.  6,  gust  suppression  by  plunging 
airfoil  is  described.  Conclusions  are  drawn  in  Sec  7. 

II.  Governing  Equations 

The  governing  equations  are  the  2-D  compressible  Navier-Stokes  equations.  To  keep  the  formulation 
applicable  to  unsteady  aerodynamics  of  flapping  airfoils,  the  Navier-Stokes  equations  are  cast  in  moving 
grid30. 
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(1) 


dU  dE  dF 
dt  dx  dy 


=  0 


On  introducing  a  transformation  (x,y, t)->(<p(x,y, t), y/(x,y, t), t),  the  governing  Navier-Stokes  equations 
become 


dU  dtp  dU  dy/  dU  dtp  dE  dy/  dE  dtp  dF  dy/  dF  v 

- +— - +  — - - +— —  +  — - +— —  +  — - =  S(t) 

dr  dt  dtp  dt  dy/  dx  dtp  dx  dy/  dy  dtp  dy  dy/ 


(2) 


where  U  denotes  the  solution  vector,  E  and  F  are  the  fluxes  (includes  both  the  inviscid  and  viscous  flux 
components)  and  ( tp ,  y/)  denotes  the  potential  function  and  streamfunction  respectively.  S  (t)  refers  to 
acoustic  or  vortical  sources  of  sound  that  are  presented  in  the  flow. 

The  solution  vector  and  the  flux  components  are  given  by 


U  = 


P 

pu 

pv 


(3) 


E  - 


pu 

pu2  +P-T „ 
P“V  -  Tr, 

(E,  +  p)u  -  UT^  -  vz^ 


+  qx 


(4) 
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(5) 


pv 

P“V  ~  *xy 

PV2+P-ryy 

(E,  +  p)u-urxy-VTyy+qy 


And 


E,  =  p(e  +  0.5(u2  +  v2)) 


(6) 


The  shear  stresses  and  the  heat  flux  vector  are  given  by 


2  du  dv 

Txx  = - (2 - ) 

“  3  Re  dy 


2  dv  du 

Tyy  3  Re '  dy  dx 


1  du  dv 

= (2  —  +  — ) 

Re  dy  ch: 


= 


-1 


dT 


(y  -  \)M  RePr  dx 


Vv  = 


-1 


dT 


(y  - 1  )M  Re  Pr  dy 


(7) 

(8) 

(9) 


(10) 


(11) 


where  Re,  M,  Pr  are  the  Reynolds  number,  Mach  number  and  Prandtl  number.  The  equations  of  state 
for  ideal  gas  become 

p  =  (y-\)pe  (12) 
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(13) 


T  =  yM^p_ 

P 

The  grid  metrics  are  calculated  using  the  following  equations 


dx  dy 
dtp  _  dx 

1/ 

dy  dy 

dY  _  jdy 

dx  d<f> 

dy  _  jdx 

dy  dtp 

dtp  _  dx  dtp  dy  dtp 

dt  dr  dx  dr  dy 

dy  _  dx  dp  dy  dtp 
dt  dr  dx  dr  dy 


(15) 

(16) 

(17) 

(18) 

(19) 


where  J  is  the  Jacobian  of  the  transformation  which  is  given  by  J  = - .The  grid 

(**>v  ~x¥y<p) 

velocities  (xt,  yT)  are  computed  from  analytical  expression  for  the  rigid  body  motion  of  the  grid  around 
the  moving  body.  The  grid  metrics  are  computed  using  the  same  high-order  spatial  discretization  scheme 
(see  Section  IV)  used  for  the  computation  of  the  fluxes. 
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III.  Grid  Generation 


Streamfunction  and  potential  function  are  conjugate  harmonic  functions  and  are  orthogonal  to  each 
other.  This  property  can  be  utilized  for  the  generation  of  orthogonal  numerical  conformal  grids.  This  idea 
has  been  employed  for  2D  grid  generation  and  in  some  cases  extended  for  3D  problems.  For  simple 
geometries,  analytical  solutions  can  be  constructed  for  streamfunction  and  potential  function  using 
elementary  solutions  of  Laplace  equation.  For  such  geometries,  the  exact  solution  is  a  system  of  non¬ 
linear  equations  which  can  be  inverted  using  Newton's  method  to  obtain  the  values  of  the  coordinates  (x, 
y)  when  the  values  of  streamfunction  and  potential  function  are  specified.  This  method  is  referred  in  the 
current  study  as  analytical  grid  generation.  For  problems,  where  the  exact  solution  is  not  available,  the 
numerical  values  of  stream  function  and  potential  function  on  the  solid  body  can  be  computed  using  panel 
methods.  The  numerical  values  of  streamfunction  and  potential  function  on  the  solid  body  are  used  to 
distribute  streamlines  and  potential  lines  in  the  numerical  domain  and  an  approximate  system  of  non¬ 
linear  equations  for  streamfunction  and  potential  function  can  be  inverted  to  obtain  the  coordinates 
(x,y).This  method  is  referred  as  the  numerical  grid  generating  method.  For  problems  involving  multiple 
bodies  in  the  domain,  it  is  possible  to  generate  the  grid  by  applying  the  panel  method  for  the  system  of 
bodies  directly  or  by  applying  the  panel  method  separately  to  each  body  and  joining  the  separate  grids 
using  stitching. 

Potential  lines  become  coarse  close  to  the  stagnation  points.  For  accurate  numerical  representation  of 
the  problem,  this  issue  needs  to  be  resolved.  In  the  current  study,  a  hyperbolic  sine  function  is  used  to 
cluster  the  potential  lines  and  it  is  found  that  the  grid  coarse  elements  close  to  the  stagnation  points  can  be 
removed.  A  hyperbolic  tangent  function  is  employed  to  cluster  the  streamlines  for  resolving  the  boundary 
layer. 

A.  Analytical  Method 
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The  analytical  grid  generating  system  is  the  faster  of  the  two  above-listed  methods  as  it  involves  just  a 
single  step  inversion  of  the  non  linear  system  of  equations.  The  deviation  from  orthogonality  for  grids 
generated  using  the  analytical  method  is  negligible.  The  analytical  grid  generation  method  is  implemented 


as  follows: 


•  Construct  the  exact  analytical  solution  for  streamfunction  and  potential  function. 

•  Generate  the  numerical  domain  between  certain  maximum  and  minimum  values  of 
streamfunction  and  potential  function. 

•  Find  (x,  y)  coordinates  for  any  given  grid  point  (yK  vj).  To  achieve  it,  invert  the  non-linear  system 
of  equations  using  the  Newton's  method  within  a  specified  tolerance.  For  the  left  comer  point  a 
random  initial  guess  is  used.  For  rest  of  the  grid  points,  the  previous  grid  point  was  used  as  an 
initial  guess. 

•  For  boundary  layer  resolution,  streamlines  are  clustered  using  the  following 

formula:  5  = - 1= ,  where  N  is  the  number  of  points  across  the  boundary  layer  and  i _ is 

Aw  Re  VRe 

the  scale  of  thickness  of  the  boundary  layer. 

An  example  of  the  grid  generated  using  the  above  algorithm  is  considered  for  the  case  of  uniform  flow 
past  a  single  cylinder  (see  Fig.  la).  For  this  case,  the  analytical  expressions  for  streamfunction  and 
potential  function  are  given  by 


(20a) 
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v  =y(  i- 


R 


x 2  +  y 2 


(20b) 


As  can  be  seen  in  Fig.  la,  the  uniform  grid  obtained  is  coarse  close  to  the  stagnation  point.  The  way  of 
avoiding  this  is  discussed  in  the  section  D. 

B.  Numerical  Grid  Generating  System 


The  numerical  grid  generating  system  consists  of  two  parts:  i)  apply  a  panel  (boundary  element 
method)  solver  to  find  the  velocity  the  potential,  and  ii)  invert  an  approximate  numerical  expression  for 
streamfunction  and  potential  function  to  find  (xt,yj)  coordinates  for  any  given  grid  point  {<p,,v ,  ).The 

orthogonal  property  of  the  numerical  method  is  comparable  to  the  hyperbolic  grid  generating  method  with 
the  advantage  of  being  able  to  specify  the  outer  boundary.  The  following  algorithm  is  used  for  the 
numerical  grid  generation: 


•  Construct  a  panel  method  solution  of  the  potential  function  on  the  solid  boundary.  In  the 
current  study,  discrete  potential  vortices31  are  used  as  an  elementary  solution  in  panel 
methods.  The  panel  method  solves  for  the  discrete  potential  vortex  strengths  o,  at  each 
boundary  cell  (panel). 

•  Since  the  surface  of  the  body  is  a  streamline,  the  ideal  value  of  streamfunction  will  be 
constant  on  the  surface  of  the  body.  The  value  of  streamfunction  on  the  surface  of  the  body 
computed  using  the  panel  methods  will  have  some  numerical  error  from  the  computations. 
Hence  the  numerical  value  of  the  streamline  will  not  be  the  same  at  all  the  points  on  the 
surface  of  the  body.  To  account  for  this  error,  the  values  of  the  streamline  on  the  body  are 
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denoted  as  Cu  and  Q  where  Cu  is  the  maximum  value  of  the  streamline  on  the  body  and  Q  is 
the  minimum  value  of  the  streamline  on  the  body. 

•  Next,  the  first  streamline  above  and  below  the  body  are  constructed  as  \|/,=CiCu  and  y.i=C2C|, 
where  Ct,  C2  are  used  to  ensure  that  the  streamlines  does  not  cross  into  the  body.  The  values 
of  C|  and  C2are  set  to  1.02  and  0.98  respectively.  Further,  to  ensure  that  the  required  number 
of  points  is  clustered  in  the  boundary  layer,  the  following  condition  is  satisfied  in  the 

numerical  algorithm:^  =|  V'i  ~¥-\  I  ,  where  8  = - \=  .The  number  of  discrete  vortices 

Aw  Re 

(panels)  is  increased  on  the  surface  of  the  solid  body  and  the  panel  method  is  repeated  until 
the  above  condition  is  satisfied. 

•  Generate  the  numerical  domain  between  certain  maximum  and  minimum  values  of 
streamfunction  and  potential  function. 

•  The  following  system  of  equations  is  inverted  to  obtain  the  values  of  (x,y) 


<!>(x,y)  =  x  - 1  tan'1  (-  y' ) 
2k  x-x, 


(21a) 


(21b) 


An  example  of  the  application  of  the  above  method  is  shown  for  the  case  of  a  cambered  airfoil  (Fig.2). 
It  can  be  seen  that  the  grid  lines  are  orthogonal  to  each  other.  The  above  grid  generating  algorithm  has  the 
ability  to  produce  high  quality  orthogonal  grids  on  different  kinds  of  geometries  (Fig.  1-4). 
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Elliptic  grid  generating  system  can  generate  smooth  grids  but  the  implementing  of  boundary  layer 
clustering  and  generating  orthogonal  grids  require  some  effort.  The  hyperbolic  grid  generation  can 
generate  orthogonal  grids  but  the  outer  boundary  cannot  be  specified  and  the  grid  can  deteriorate  rapidly 
close  to  sharp  edges  unless  grid  lines  are  clustered  or  artificial  dissipation  is  added.  The  current  grid 
generation  method  can  be  viewed  as  a  hybrid  method  which  is  able  to  achieve  the  orthogonality  and  grid 
clustering  with  the  ability  to  specify  the  outer  boundary  at  the  same  time.  The  implementation  of  grid 
clustering  will  be  discussed  in  Sec. III. D. 

C.  Grid  Stitching 

For  problems,  involving  multiple  bodies  like  two  wings  in  tandem  or  a  combination  of  steady  and 
unsteady  wings,  the  panel  method  can  be  used  to  generate  a  single-block  structured  grid.  An  example  is 
shown  in  Fig.3  for  the  generation  of  grids  around  two  cylinders  placed  in  tandem.  The  panel  method  grid 
generator  can  be  used  for  inviscid  studies  but  is  not  a  convenient  grid  generator  for  viscous  flows  due  to 
the  lack  of  local  grid  control  around  each  body.  This  becomes  extremely  important  for  studying  flows 
with  boundary  layer  resolution.  To  have  more  control  on  the  grid  about  individual  bodies  stitching 
method  is  employed.  In  this  method,  the  grids  around  each  body  are  generated  using  either  analytical  or 
numerical  grid  generating  system.  The  local  grids  are  then  stitched  together  at  an  interface  (stitching 
boundary) 32.  In  the  current  study,  an  elliptic  grid  generating  method  is  used  for  stitching  the  two  parts  of 
the  grid  together  to  keep  the  procedure  simple32, 33.  The  following  steps  are  used  for  generating  the 
patched  grids: 

•  Generate  the  local  streamfunction  grid  around  each  body. 
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The  region  close  to  the  patched  area  might  have  grid  lines  crossing  each  other.  To  remove 
those  grid  lines  and  make  the  grid  smooth,  the  grid  lines,  which  do  not  intersect  at  both  the 
ends,  are  selected  and  fixed. 


•  For  the  top  and  bottom  lines,  the  corners  from  the  previous  grid  lines  are  used  to  distribute 
points  on  the  comers  using  linear  interpolation. 

•  Once  the  four  comer  grid  lines  are  fixed,  an  elliptic  grid  generation  method  is  used  for 
smoothing  the  patched  region. 

The  above  steps  ensure  that  the  grid  stitching  is  local  and  does  not  change  the  global  grid.  The  size  of 
the  patching  domain  is  problem-dependent  and  needs  some  trial  and  error.  This  approach  has  been  applied 
in  the  current  study  to  generate  a  single  non-uniform  structured  grid  around  two  cylinders  in  tandem.  The 
grid  is  generated  locally  around  each  cylinder  and  is  then  smoothed  using  an  elliptic  grid  generator.  Grid 
generated  by  panel  method  around  two  cylinders  in  tandem  is  shown  in  Fig.4a.  The  grid  generated  by 
stitching  is  shown  in  Fig.4b.  The  grids  are  similar  to  each  other  if  no  clustering  is  implemented. 

D.  Grid  Clustering 

Vinokur34  found  that  the  hyperbolic  sine  functions  are  the  best  for  concentrating  points  in  the  middle 
of  domain  and  hyperbolic  tangent  function  for  clustering  at  the  ends  of  domain.  As  can  be  seen  from 
Fig.  la,  the  grid  becomes  coarse  close  to  the  stagnation  point  of  the  cylinder.  Hence  to  remove  the  coarse 
grid  cell  in  the  <p  direction,  which  is  near  stagnation  points  in  the  middle  of  the  domain,  a  hyperbolic  sine 
function  is  employed35. 


(22a) 
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where  xc  is  the  location  in  physical  domain  around  which  the  grid  clustering  is  required  (stagnation 
points  have*c  =±0.5),  rj  is  the  normalized  computational  coordinate  which  varies  from  0  to  1.  The 

clustering  parameter  x  varies  from  zero  to  infinity,  where  r =0  corresponds  to  a  uniform  grid  (no 
clustering).  The  parameter  B  is  given  by 


B=  1  l  +  (er-l  ){xjh) 
It  l  +  (e-r-l  )(xc/h) 


(22b) 


where  h  is  the  length  of  physical  domain  in  either  positive  or  negative  direction. 

For  problems  involving  viscous  flows,  resolution  of  boundary  layers  become  very  important.  To 
cluster  grid  in  the  v) /  direction  in  boundary  layer  the  following  approach  is  used13: 


*  =  Pt7  +  (1- />)(!- 


tanh[g(l-?7)] 

tanh(0 


) 


(23a) 


x  =  x,  +s(xe-x() 


(23b) 


The  factors  P  and  Q  control  the  clustering  and  x ,  and  xe  are  the  limits  of  physical  domain.  The 
comparison  between  grid  generated  with  the  application  of  hyperbolic  sine  function  (22)  is  shown  in 
Fig.  lb  with  a  clustering  parameter  r  =5.  The  clustering  of  points  in  the  boundary  layer  for  the  SD  7003 
airfoil  is  shown  in  Fig.2b  with  P=0.I  and  Q=4  in  Eq.  (23),  and  x=2.5  in  Eq.  (22). 

E.  Orthogonality  in  2D 

The  orthogonality  property  of  the  grid  is  measured  by  computing  the  angle  between  grid  lines,  which 

u 

is  given  by 
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(24) 


0  =  cos-1  (  .  ==  ) 

"V  &II&22 

where  g  is  the  metric  tensor  which  is  calculated  by 

l,dxk.,dxk  . 

gij  =  (25) 

dq,  d<Jj 


where  (x,,x2)  =  (x,.y)  and(<7,,<72)  =  (<f>,i//) .  For  an  orthogonal  grid,  0  should  be  equal  to  90  degrees. 
The  maximum  measured  deviation  from  orthogonality  was  around  4-5  degrees  for  both  grids  about  the 
cylinder  and  airfoil.  The  deviation  was  observed  only  at  two  points  to  the  front  and  back  stagnation  point 
of  the  cylinder  and  at  the  leading  edge  of  the  airfoil.  At  the  rest  of  the  grid  points  the  deviation  was  less 
than  0.5  degrees.  The  orthogonal  property  of  the  current  grid  generating  method  was  compared  with  the 
other  grid  generating  methods  available  in  literature.  Comparison  between  the  different  methods  is  shown 
in  Table.  1.  The  deviation  from  orthogonality  for  the  proposed  SFC-based  non-uniform  grid  appears  to  be 
comparable  to  other  listed  methods. 


Table  1.  Comparison  of  orthogonal  property  of  different  orthogonal  grid  generating  system 


Method 

Maximum  Departure  (in  degrees) 

Geometry 

Steger  and  Chausse36 

16 

Turbine  blade 

Duraiswamy  and  Prosperetti37 

1.11 

NACA  0015 

Eca38 

4.41 

NACA  0015 

Nair  and  Sengupta39 

0.5 

NACA  0015 

Present 

4-5 

SD  7003/NACA  0012/NACA 
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IV.  Numerical  Algorithm 


0015 


The  governing  equations  are  discretized  using  high-order  finite-difference  schemes.  This  is  done  to 
reduce  significantly  the  number  of  points  required  per  wavelength  compare  to  second-order  schemes  for 
the  resolution  of  the  acoustic  waves  and  vortices24. 

A.  Spatial  Differences 

A  non-uniform  compact  scheme  was  proposed  in  Ref.  25.  The  scheme  becomes  identical  to  the 
uniform  compact  scheme  proposed  by  Ref.24  .For  some  function /its  derivative /  is  computed  as 

+  fi  +  Pi  fi+\  =  ^J^ilfi+I  (26) 

/=- 2 

The  coefficients  on  the  left-hand  side  and  right-hand  side  can  be  obtained  by  matching  the  Taylor’s  series 
expansion.  Close  to  the  boundaries  at  i=l,  2  and  i=  N-l,  N,  the  scheme  is  replaced  with  one-sided 
schemes  and  the  coefficients  are  obtained  accordingly.  More  details  can  be  found  in  Refs.  14,  15.  The 
value  of  or,  and  /?,  can  be  adjusted  at  each  point  but  are  kept  constant  at  1/3  at  each  point  and  the 
scheme  becomes  central-differenced.  On  a  non-uniform  grid,  the  fourth-order  accuracy  is  achieved. 


B.  Numerical  filtering 
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Since  the  spatial  differencing  scheme  is  central-difference,  numerical  instabilities  arising  due  to 
unresolved  scales  and  grid  non-uniformities  can  grow.  Hence  the  solution  is  filtered  according  to  the 

procedure  suggested  in  Ref.  16.  If /is  the  unfiltered  solution  and  f  is  the  filtered  solution  then 

A  5 

affi- 1  +  fi  +affi+\  =  ^ ailfi+l  (27) 

/=- 5 

Here  af  is  a  free  parameter  with  values  between  -0.5  and  0.5  based  on  analysis  using  transfer 

function.  The  higher  the  value,  the  less  dissipative  is  the  filter.  In  the  current  study  the  value  is  set  equal 
to  0.45. Since  the  stencil  is  wide,  one-sided  differencing  formulas  are  needed  at  points  (1 ...  5)  and  (N-4... 
N).  Higher-order  one-sided26  formulas  can  be  used  to  obtain  the  filtered  values  close  to  the  boundaries. 
Ref  26  has  a  detailed  discussion  on  the  derivation  and  the  order  of  the  filters.  The  numerical  filter 
employed  in  the  current  study  is  the  tenth-order  accurate. 

C.  Time  Marching 

Low-storage  Runga-Kutta  proposed  by  Ref.27  is  used  for  marching  the  solution  in  time.  The  method 
is  explicit  and  third-order  accurate.  The  time  marching  for  a  system  given  by  df/dt=u  is  computed  from 
time  step  m  to  time  step  m+1  as  follows 

fn+'=fn+bn+'Hn  (28) 

Hn  =  a»Hn-\  +  dtU  ,  (29) 

where  n  refers  to  the  stage  number.  The  value  of / at  the  final  stage  of  the  scheme  becomes  the  value  at 
the  timestep  m+1.  The  values  of  the  coefficients  a„  and  bn+/  have  been  given  in  Ref.27. 
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D.  Boundary  Conditions 


Different  kinds  of  boundary  conditions  have  to  be  specified  for  the  Euler  and  Navier  Stokes  equations. 
For  Euler  equations  at  the  wall,  ghost  point  boundary  condition  of  Tam  and  Dong28  in  employed.  No-slip 
boundary  conditions  are  used  for  the  stationary  Navier-Stokes  equations  and  the  surface  velocity 
components  are  set  to  the  grid  speed  for  the  moving  grid  case.  At  the  artificial  domain  boundary  energy 
transfer  and  annihilation  (ETA)29  is  employed  with  extrapolation  of  the  values  to  domain  boundaries. 


V.  Numerical  Results 


A.  Euler  Solver 

As  a  first  step,  the  code  is  validated  for  non-linear  and  linearized  Euler  equations  (Eq.2,  without  the 
viscous  terms)  for  curvilinear  geometries  using  the  SFC.  A  problem  from  the  2nd  CAA  workshop40  is 
considered.  The  configuration  describes  the  scattering  of  an  initial  acoustic  pulse  from  a  circular  cylinder 
(See  Fig.5a).The  initial  acoustic  pulse  is  given  by 

-ln(2X(x-*c)2+(.y->02) 

p  =  pn  +  O.Ole  004  (30) 


The  center  of  the  pulse  (xc  ,yj  is  located  at  (4,0).  The  domain  is  kept  uniform  in  the  region  [-8*8]*[- 
8*8]  and  coarsened  beyond  that  with  a  stretching  factor  of  1.2.  The  results  were  computed  on  300*300 
and  450*450  grids  and  the  results  were  found  to  converge.  The  pressure  contour  for  pulse  reflection  from 
the  cylinder  and  propagation  into  region  of  grid  stretching  are  shown  in  Figs.6a,  b.  As  can  be  seen  in 
Fig.6b,  the  grid  stretching  and  filtering  absorbs  the  outgoing  energy  without  producing  non-physical 
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reflections  back  into  the  domain.  The  results  are  compared  on  different  grids  at  PI  (0,  0.6)  (Fig.5b).  The 
results  obtained  on  the  300*300  and  250*250  grid  agreed  with  each  other.  For  the  200*200  grid,  two  sets 
of  computations  were  performed.  It  can  be  clearly  seen  from  Fig.5b,  the  computations  with  clustering  of 
grid  points  provides  a  closer  match  with  the  results  obtained  on  the  two  previous  grids.  For  the  case  with 
no  grid  clustering  because  of  the  coarse  grid  cells  near  the  stagnation  point,  the  number  of  grid  points  per 
wavelength  is  reduced.  This  causes  smearing  of  the  acoustic  pulse.  The  clustering  increases  the  number  of 
points  close  to  the  stagnation  region  and  ensures  that  sufficient  numbers  of  points  are  present  to  resolve 
the  waves. 

Two  different  kinds  of  grid  generation  were  mentioned  in  Section  1II.C  for  problems  involving 
multiple  bodies.  To  test  the  grid  generation  algorithm  for  multiple  bodies  using  both  the  grid  generation 
methods,  sound  scattering  from  two  cylinders  placed  in  tandem  is  considered  (Fig.4).  An  initial  pressure 
pulse  is  similar  to  the  single  cylinder  case(Eq.30)  and  is  located  at  the  origin.  The  aeroacoustic 
computations  were  performed  with  the  left  cylinder  of  diameter  1  located  at  (-4,  0)  and  the  right  cylinder 
of  diameter  0.5  located  at  (4,  0).The  domain  of  computation  is  [-8,  8]*[-5,  5].  The  acoustic  pressure  was 
monitored  at  the  point  (0,  2).  As  shown  in  Fig.7a  excellent  agreement  is  obtained  between  both  methods 
of  grid  generation. 

To  test  the  long  term  stability  of  the  solver  and  the  boundary  conditions,  the  scattering  of  periodic 
pulse  placed  from  two  cylinders  in  tandem  is  considered41.  The  aeroacoustic  computations  were 
performed  with  the  same  two  cylinder  tandem  configuration  but  only  on  the  stitched  grid.  A  periodic 
source  given  by 

5  =  0. Ole  004  sin(ry/)  (31) 


was  located  at  (0,  0).The  value  of  w  was  taken  to  be  8jt.  The  periodic  source  was  added  to  the  right-hand 
side  of  the  energy  equation  and  continuity  equation.  The  initial  values  of  the  acoustic  pressure  and 
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velocities  were  assumed  to  be  zero.  The  region  close  to  the  cylinder  surface  was  stretched  using  a 
clustering  parameter  r  =1  in  Eq.  (22).  The  degree  of  stretching  was  kept  limited  to  ensure  at  least  8  points 
per  wavelength  in  the  far  field.  The  computations  were  performed  on  the  720*720  grid.  The  pressure 
history  at  a  probe  point  is  shown  in  Fig.7b  and  it  can  be  seen  that  once  the  transient  waves  leave  the 
domain,  the  solution  has  become  periodic.  The  rms  value  of  pressure  was  compared  on  the  surface  of  the 
left  (Fig.8a)  and  right  cylinder  (Fig.8b)  with  the  analytical  solution41.  Very  good  match  is  obtained  for 
both  the  cases  and  demonstrates  the  accuracy  of  computation  for  complex  geometries. 

B.  Navier-Stokes  Solver  in  Fixed  Coordinates 

To  test  the  Navier-Stokes  solver  on  stationary  grids,  the  problem  of  uniform  flow  past  a  circular 
cylinder  is  considered  first  (see  Figs.  9-10).  This  test  case  was  chosen  since  abundant  experimental  and 
computational  data  are  available  for  comparison  of  shed  vortex  pattern  in  the  literature.  It  has  been  found 
that  the  vortex  shedding  appears  for  Re  around  40-50  as  observed  in32  and  the  references  mentioned 
therein.  Below  the  value  of  Re  equal  to  40,  a  pair  of  symmetric  vortices  is  formed  in  the  wake  of  the 
cylinder.  The  size  of  the  vortices  increases  with  Re  until  a  value  of  Re=50  is  reached,  when  the  periodic 
vortex  shedding  begins. 

All  the  simulations  were  performed  with  the  grid  sizes  of  100*100,  150*150  and  200*200.The  results 
obtained  with  the  150*150  and  200*200  grids  are  coincided.  To  resolve  the  viscous  boundary  layer,  grid 
stretching  was  employed  by  adjusting  the  parameters  P  and  Q  to  provide  at  least  15-20  points  in  the 
boundary  layer.  It  was  found  that  for  the  cases  considered  no  grid  clustering  was  needed  in  the  (p 
direction. 

The  computations  were  performed  at  Re=20,  40,100  and  150.  Table  2  presents  a  comparison  of  the 
current  simulations  and  the  available  data.  Good  agreement  is  found  at  these  Reynolds  numbers.  For  the 
Re=20  and  Re=40  (Figs.9b,  c),  steady  symmetric  vortices  are  formed  behind  the  cylinder.  The  size  of  the 
vortices  was  compared  with  the  available  data  from  literature32.  For  the  cases  of  Re=100  and  150,  once 
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the  transient  flow  leaves  the  domain,  periodic  vortex  shedding  appears  behind  the  cylinder.  Table  2 
compares  the  Strouhal  number  (St)  for  vortex  shedding  obtained  in  the  current  study  with  results  that  has 
been  published  in  the  literature.  Again  the  good  agreement  is  found  in  terms  of  the  Strouhal  number. 


Reynolds’s  Number 

Current  Simulations 

Literature 45 

20 

L=0.92  a=0.34 

L=0.9 1-0.93 

b=0.42 

a=0.33-0.36 

b=0.43-0.46 

40 

L=2. 13-2.28 

L=2.15 

a=0.71-0.76 

a=0.7 

b=0.59-0.6 

b=0.6 

100 

St=0.167 

St=0. 1 6-0. 1 7 

150 

St=0. 181 

St=0. 18-0. 187 

Table  2.  Comparison  of  the  vortex  patterns  between  current  simulation  and  data  available  in 
the  literature 


To  test  the  Navier-Stokes  solver  in  the  range  of  moderate  Reynolds  numbers  (10000-60000),  uniform 
flow  past  a  stationary  airfoil  (SD  7003)  was  studied  at  a  Re=20000  over  a  range  of  angle  of  attacks  (0-12 
degrees).  It  was  found  that  using  the  very  mild  value  of  grid  clustering  parameter  r =0.1  the  computed 
natural  frequency  of  vortex  shedding  (f=11.5)  was  close  to  the  measured  frequency  of  vortex  shedding 
(f=12-13). 


C.  Navier-Stokes  Solver  in  moving  coordinates 
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U  c 

The  aerodynamics  of  a  plunging  airfoil  with  Re  =  — —  =  20000  was  studied  as  a  test  case  for  the 

v 

Navier-Stokes  solver  on  the  moving  grid.  The  simulations  were  performed  for  NACA  0012  airfoil  and  SD 
7003  airfoil. 

The  plunging  motion  of  an  airfoil  is  given  by 

y(t)  =  hcos(2kt)  (32) 

coc 

where  h  is  the  plunging  amplitude  normalized  by  the  chord  length  and  k= -  is  the  reduced 

2(7 

oo 

frequency.  The  NACA  0012  airfoil  case  was  run  for  the  Strouhal  number  St=kh=0.3.  The  computations 
were  performed  on  four  grid  levels  (G1  (240*80),  G2  (300*80),  G3  (360*80)  and  G4  (400*800)).  The 
grids  G3  and  G4  were  generated  with  grid  clustering  parameter  r  =  1  in  Eq.  (22).  The  G2  was  generated 
with  the  grid  clustering  parameter  r  =2.0  in  Eq.  (22).  After  an  initial  transient  period  of  2  cycles,  the 
results  become  periodic.  The  results  on  G3  and  G4  coincide  while  the  results  on  G1  showed  smearing. 
Peak  coefficient  of  lift  obtained  using  different  grids  (G2  and  G3)  and  the  results  of  Young43  are  shown  in 
Fig.  1 1 .  The  obtained  computational  results  correspond  to  those  of  Young43  for  the  wide  range  of  it. 

The  aerodynamics  of  a  plunging  SD  7003  airfoil  were  performed  for  two  sets  of  data  (i)  h=0.025,  k= 
7.85  (St=0.196)  and  (ii)  h=0.1,  k  =  5.91  (St=0.591).  The  former  case  is  designated  as  a  low  Strouhal 
number  case  and  the  latter  case  is  a  high  Strouhal  number  case.  The  grid  speedy  is  calculated  from 

Eq.32.  The  computations  were  performed  using  G2  and  G2  with  T  =3.5.  It  was  found  that  for  the  low  St 
case  the  computations  ran  on  both  grids  without  producing  any  unphysical  reflections.  For  the  high  St 
case,  the  computations  with  mild  clustering  produced  large  amounts  of  unphysical  oscillations  close  to  the 
leading  edge  which  contaminated  the  solution.  To  remove  the  unphysical  oscillations,  the  clustering 
parameter  was  changed  to  3.5  and  the  unphysical  oscillations  were  avoided.  For  high  St,  strong  leading 
edge  vortex  is  observed  as  shown  in  Fig.  12. 
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The  flow  field  was  compared  with  the  flow  visualization  of  Ol44  for  both  the  cases  by  plotting  the 
hyperbolic  tangent  of  vorticity.  Similar  vorticity  pattern  to  the  experiments  were  observed  during  the  top 
stroke  and  bottom  stroke  (Fig.  13)  showing  the  accuracy  of  the  computations. 

VI.  Gust  Suppression  by  Flapping  Airfoils 

In  the  current  study,  the  effect  of  the  wind  gust  on  the  aerodynamic  performance  of  a  2D  plunging 
airfoil  is  examined  numerically  using  a  Navier-Stokes  solver  in  moving  coordinates  described  above.  The 
gust  is  included  as  a  source  term45  in  the  right-hand  side  of  the  momentum  equations  to  avoid  dealing 
with  the  numerical  problems  to  propagate  gust  at  the  boundaries  of  the  computational  domain.  The  flight 
of  Micro  Air  Vehicles  (MAV)  experiences  high  unsteady  flow  conditions  due  to  the  gusty  winds  in  the 
atmosphere.  Birds  and  insects,  which  operate  under  similar  conditions,  use  the  flapping  motion  of  their 
wings  to  generate  thrust  and  lift.  To  improve  the  design  of  the  MAV  to  suit  the  operating  conditions,  a 
number  of  studies  have  been  performed  to  study  the  aerodynamics  of  insects  and  simplified  models  of 
airfoils.  However,  the  effect  of  the  wind  gust,  which  is  a  very  important  factor  in  the  MAV  design,  has 
been  neglected  in  most  of  these  studies.  In  the  current  study,  the  effect  of  the  wind  gust  on  performance  of 
MAV  and  the  ability  to  suppress  the  effect  of  wind  gust  using  the  plunging  motion  of  the  airfoil  is 
examined. 

The  gust  is  included  as  a  source  term  in  the  right-hand  side  of  the  momentum  equations  to  avoid 
dealing  with  the  problems  at  the  boundaries  of  the  computational  domain: 

(  °  ) 

s=  C|  cos(kt) 

C2  cos(Af) 

<  0  j 

wc 

where  k  is  the  reduced  frequency  given  by  k  = - .  The  coefficients  C, ,  C2  are  given  by 

2£/m 

C,  =  3a(l  +  cos (b(x  -  xc  )))(tanh(3(>'  +  yc  ))2  -  tanh(3(>>  -  yc ))2 )  (33) 
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C2  =  ab  sin  (b(x  -  xc  )))(tanh(3(^  +  yc ))  -  tanh(3(>>  -  yc ))) 


(34) 


width,  xc,yc  represent  the  center  of  the  vortical  gust. 


The  amplitude  of  plunging  (h)  is  computed  from  the  formula 


h  =  e  max(C,  ,C2), 


(35) 


where  £  is  a  parameter,  which  is  used  to  adjust  the  plunging  amplitude  to  match  the  lift  force 
generated  by  the  vortical  gust  with  the  lift  force  generated  by  the  plunging  motion  of  airfoil. 

Computations  were  performed  for  three  cases:  (plunging  frequency  and  the  vortical  gust  frequency 
were  equal) 

o  Stationary  airfoil  with  vortical  gust  (Case  1) 
o  Plunging  airfoil  with  no  vortical  gust  (Case  2) 

o  Plunging  airfoil  used  to  suppress  the  oscillatory  lift  force  created  by  the  vortical  gust.(Case  3) 
The  comparisons  are  shown  in  Fig.  1 4  fora  =  0.8, b  =  5 ,h  =  0.007, £  =  3.93  .  As  the  gust  propagates  with 
the  mean  flow,  the  approximate  non-dimensional  wavelength  of  the  gust  is  0.8.  As  the  non-dimensional 
chord  length  of  the  airfoil  in  the  computation  is  equal  to  unity,  more  than  one  wavelength  of  the  gust  is 
swept  across  the  airfoil  during  each  instant.  The  computations  for  case  1  and  case  2  in  Fig.  1  shows  that 
the  case  1  and  case  2  are  approximately  out  of  phase.  For  the  case  3,  it  is  found  that  lift  suppression  has 
been  achieved  but  the  suppression  is  not  a  simple  superposition  of  the  lift  force  generated  due  to  the  gust 


and  the  plunging  motion.  For  more  exact  cancellation  of  the  gust  lift,  the  pure  plunging  should  be 
combined  with  the  pitching  motion  of  airfoil  that  will  be  considered  in  our  future  research. 

VII.  Conclusions 

Streamline  as  a  Coordinate  method  (SFC)  was  used  in  the  past  for  steady  aerodynamics  with  low- 
order  of  approximation  of  governing  partial  differential  equations.  The  advantages  of  the  SFC  are  the 
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ability  of  grid  to  adhere  to  the  curvilinear  airfoil  surface,  to  resolve  the  boundary  layer,  and  to  move  with 
the  pitching  and  plunging  airfoil.  The  use  of  SFC  for  aerodynamics  of  flapping  airfoils  was  hindered  by 
low  accuracy  of  the  SFC  method  near  stagnation  points. 

In  the  current  study  the  SFC  method  using  analytical  and  numerical  solutions  of  potential  function  and 
patching  of  sub-domains  was  combined  with  the  high-order  spatial  and  temporal  discretizations  in  a 
single-grid  approach.  The  far-field  SFC  grid  approaches  Cartesian  grid  so  the  ETA  method  was  found  to 
absorb  the  outgoing  numerical  error.  The  presented  results  show  the  orthogonality  property  of  the  grid 
generated.  This  paves  the  way  to  adopt  the  economic  single-grid  approach  for  complex  geometry  when 
high-order  schemes  are  needed.  However,  it  was  observed  that  a  uniform  distribution  of  SFC  grid  points 
in  the  computational  domain  produces  coarse  grid  cells  close  to  the  stagnation  point  that  may  harm  the 
effort  of  utilization  of  high-accuracy  high-order  schemes.  To  circumvent  this  issue,  a  sine  hyperbolic 
stretching  function  was  used  for  grid  refinement  in  the  (p  direction  in  the  current  study. 

The  numerical  solution  of  governing  Navier-Stokes  in  non-inertial  coordinates  have  shown  that  the 
non-clustered  grid  in  the  9  direction  in  the  vicinity  of  stagnation  points  can  be  used  if  the  modeled  flow  is 
steady,  that  is,  the  only  velocity  scale  is  present  and  this  velocity  tends  to  zero  near  the  stagnation  point.  If 
another  non-zero  velocity  scale  is  present  such  as  the  speed  of  flapping  airfoil,  the  grid  clustering  is 
needed  in  the  9  direction  in  the  vicinity  of  stagnation  points. 

The  scattering  of  sound  from  single  and  multiple  cylinders  was  taken  as  representative  cases  for  which 
analytical  or  well-established  numerical  solutions  are  available  for  inviscid  Euler  equations.  For  the  case 
of  a  single  cylinder  and  multiple  cylinders,  the  computed  acoustic  pressure  matched  well  with  the  exact 
solution  when  appropriate  clustering  of  grid  points  is  implemented.  Uniform  flow  past  a  single  cylinder 
was  taken  as  a  test  case  for  viscous  flows  at  low  Reynolds  numbers.  The  size  of  shed  vortices  for  this  case 
had  a  good  agreement  with  the  data  available  in  the  literature.  For  the  periodic  shedding  at  Re=100  and 
150,  the  Strouhal  number  matched  with  the  values  given  in  the  literature.  The  grid  refinement  in  the  9 
direction  in  the  vicinity  of  stagnation  points  was  not  needed  since  the  only  velocity  scale  is  present. 
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The  vortex  shedding  from  a  SD  7003  airfoil  in  uniform  flow  was  studied  as  a  test  case  for  moderate 
Reynolds  number  flow  (Re=20000)  and  it  was  found  that  the  computed  natural  shedding  frequency 
agreed  well  with  the  experimental  results  with  very  mild  grid  clustering  in  the  (p  direction. 

Plunging  of  airfoils  at  two  different  Strouhal  numbers  was  studied  to  find  the  necessary  degree  of  grid 
clustering  near  stagnation  points.  The  measured  coefficient  of  lift  for  the  case  of  NACA  0012  airfoil 
coincided  with  the  results  obtained  by  Young43.  The  flow  field  visualization  for  SD  7003  plunging  airfoil 
showed  similar  structures  to  the  observed  experimental  pattern  when  appropriate  clustering  of  grid  points 
in  the  (p  direction  is  implemented.  Further,  it  was  found  that  the  amount  of  grid  clustering  required  in  the 
(f>  direction  increased  for  higher  Strouhal  numbers,  that  is,  for  a  higher  plunging  speed  of  the  moving 
airfoil. 
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Figure  1.  The  SFC  grid  around  a  cylinder:  (a)  uniform  grid  and  (b)  grid  stretching  in  the  <p-direction  close  to  the  stagnation 

points.(r=4,  in  Eq.  (22)) 


Figure  2.  Grid  generated  using  panel  method  for  SD  7003  Airfoil:  Domain  with  boundary  layer  v|/-clustering  and  the  tp-clustcring 

close  to  stagnation  point 
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(a) 


(b) 


Figure  3.  Grid  generated  using  panel  method  for  SD  7003  Airfoil:  (a)  zoomed-in  view  of  leading  edge  and  (b)  zoomed-in  view  of 

trailing  edge 
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Figure  4.  Grid  generation  for  pair  of  cylinders  in  tandem:  (a)  using  panel  method  solution  of  two  cylinder  system  and  (b)  local 
panel  method  solution  for  each  cylinder  followed  by  stitching  along  the  region  of  patching.  Every  second  grid  line  is  shown  on 
either  side  starting  from  the  streamline  on  top  and  bottom  of  the  cylinder  for  both  the  cases. 
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Figure  5.  Pressure  contour  for  scattering  of  an  initial  pulse  from  cylinder  a)  t=5  and  b)  t=9. 


(a) 


(b) 


Figure  6.Results  on  the  scattering  of  acoustic  pressure  from  a  single  cylinder:  (a)  numerical  grid  with  coarsening  and  (b)  acoustic 
pressure  at  PI  obtained  by  using  a  uniform  grid  (300*300)  and  a  grid  of  200*200  with  and  without  clustering  close  to  stagnation  point 


34 

Final  report  FA9550-07-1-0314 


3/13/2008 


Acoustic  Pressure 


(a)  (b) 

Figure  7.  Sound  scattering  from  the  pair  of  cylinders:  (a)  momentarily  acoustic  pulse  computed  on  grids  obtained  by  grid 
stitching  method  and  by  the  single-domain  panel  method  and  (b)  pressure  history  at  probe  point  for  scattering  of  a  periodic 

pulse 


(a)  (b) 

Figure  8.Comparison  of  analytical  and  numerical  solutions  for  scattering  of  a  periodic  pressure  pulse  from  the  pair  of  cylinders:  (a) 
acoustic  pressure  on  the  surface  of  the  left  cylinder  and  (b)  Acoustic  pressure  on  the  surface  of  the  right  cylinder 
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Nomenclature  UHd  for  Comparison  of  Steady  flow 


Figure  9.  Steady  vortex  past  a  single  circular  cylinder:  a)  nomenclature  used  for  comparing  quantities  of  interest, 

b)  Vortex  pattern  for  Re=20,  and  c)  vortex  pattern  for  Re=40. 


Figure  10.  Shedding  of  street  of  vortices  past  a  single  cylinder:  a)  Re=I00  and  b)  Re=150 
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Figure  11.  Peak  coefficient  of  lift  for  plunging  NACA  0012  airfoil 


Figure  12.  Plunging  SD  7003  airfoil  for  St=0.591:  (a)  refined  grid  and  (b)  Isolines  of  hyperbolic  tangent  of  vorticity 
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(C)  (d) 

Fig.  13  Comparison  between  experimental  observation  and  numerical  simulation  for  plunging  SD7003,  Case  1:  (a-b)  bottom 
stroke,  (c-d  )  top  stroke,  (a,  c)-experimental  visualization,  (b,  d)-computational  isolines  of  vorticity.  Hyperbolic  tangent  of 
vorticity  obtained  by  the  numerical  simulations  is  shown  in  subfigures  (b,  d) 
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Figure  14.  Comparison  of  the  Lift  force  for  the  three  cases:  Case  1  (gust),  Case  2  (plunging),  and  Case  3  (combined). 
In  the  Case  3,  airfoil  plunging  suppresses  the  effect  of  gust. 
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